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^ , Abstract 
Ph. 

, We present an algorithm to solve BSDEs based on Wiener Chaos Expansion and Picard's 

iterations. We get a forward scheme where the conditional expectations are easily computed 
thanks to chaos decomposition formulas. We use the Malliavin derivative to compute Z. 
Concerning the error, we derive explicit bounds with respect to the number of chaos and the 
discretization time step. We also present numerical experiments. We obtain very encouraging 
\ results in terms of speed and accuracy. 

Oh 

^ ! 1 Introduction 

H— > 

a 

In this paper, we are interested in the numerical approximation of solutions (Y, Z) to backward 
stochastic differential equations (BSDEs for short in the sequel). BSDEs have been introduced by 
J.-M. Bismut in [Bis73j in the linear case, whereas the nonlinear case has been considered later 
by E. Pardoux and S. Peng in [PP90j . A BSDE is an equation of the following form 

Y t=i + J t f(s,Y s ,Z s )ds- J Z s -dB s , 0<t<T, (1) 

where B is a d-dimensional standard Brownian motion, the terminal condition £ is a real-valued 
| J^-measurable random variable where {J r t}o<t<T stands for the augmented filtration of the Brow- 

• nian motion B and the generator / is a map from [0, T] x K x M. d into K. A solution to this equation 

is a pair of processes {{Y t , Z t )}o<t<T which is required to be adapted to the filtration {Tt}o<t<T- 
We will assume the conditions of Pardoux and Peng to ensure existence and uniqueness of solu- 
tions. 

X ■ 

5_j | Our main objective in this study is the numerical approximation of the solution (Y, Z) to 

BSDE |T]) (even though there exists a large literature on this subject). The first two contributions 
to this topic are due to D. Chevance |Che97j . who considered generators independent of Z, and 
V. Bally |Bal97| . who used a random time mesh. J. Ma and J. Yong }MY99| proposed numerical 
schemes based on the link between Markovian BSDEs and semilincar partial differential equations 
(PDEs). Another approach, based on Donsker's theorem and close to |Che97j . was proposed by 
F. Coquet, V. Mackevicius and J. Memin [CMM99J in the case of a generator / independent of 
Z\ the general case was treated by Ph. Briand, B. Delyon and J. Memin in [BDMOlj . who later 
extended it to a more general framework BDM02|, including the case of a "stepwise constant 
Brownian motion". This extension led to the formulas 

Y t =E(Y t + h\T t ) + hf(t, Y t , Z t ), Zt = h- 1 ' 2 E (Y t+h (B t+h - B t ) \ F t ) 
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known as the dynamic programming algorithm. Even though the convergence was proved in the 
case of path-dependent terminal condition £, the rate of convergence was left as an open question 
in |BDM02j . This problem was solved by J. Zhang |Zha04j and B. Bouchard and N. Touzi |BT04| 
in the case of Markovian BSDE, namely in the case of a terminal condition £ = g{X<r) where X 
is the solution to a stochastic differential equation. Their result was generalized by E. Gobet and 
C. Labart [GL07] and also by E. Gobet and A. Makhlouf [GMlOj . 

From a numerical point of view, the main difficulty in solving BSDEs is to efficiently compute 
conditional expectations. Several approaches have been proposed using various tools: the Malliavin 
calculus |BT04| . regression methods }GLW051 |GLW06| and quantization technics |BP03| . 

Finally, let us mention that there exists some works dealing with the discretization of solutions 
to BSDEs in a more general framework: forward-backward SDEs |DM06j and quadratic BSDEs 



Let us now describe briefly the main points of our approach in the case of a real- valued Brownian 
motion. As already used in several quoted papers, see also [BD071 IGLIOI IBSar] . our starting point 
is the use of Picard's iterations: (Y°, Z°) = (0, 0) and for q G N, 



It is well-known that the sequence (Y q ,Z q ) converges exponentially fast towards the solution 
(Y, Z) to BSDE fT]). We write this Picard scheme in a forward way 



where D t X stands for the Malliavin derivative of the random variable X . 

In order to compute the previous conditional expectation, we use a Wiener chaos expansion of 
the random variable 



More precisely, we use the following orthogonal decomposition of the random variable F q 



where Ki denotes the Hermite polynomial of degree I, (<7i)j>i is an orthonormal basis of L 2 (0, T) 
and, if n = (rii)i>i is a sequence of integers, \n\ = J2i>i n i- (^fe)fc>i,|n|=fc is t ne sequence of 
coefficients ensuing from the decomposition of F q . Of course, from a practical point of view, we 
only keep a finite number of terms in this expansion: 

• we work with a finite number of chaos, p; 

• we choose a finite number of functions gx, . . . , g^. 

This leads to the following approximation with n = (m, . . . , tin) 



One of the key point in using such a decomposition is that, for choices of simple functions g\, . . . , 
gN, there exist explicit formulas for both 



jRicllj . 








E (F q | T t ) and Z% 



■9+1 



D t E (F q | F t ) 



(2) 
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this plays a crucial role in our algorithm. Using these formulas and starting from M trajectories 
of the underlying Brownian motion we are able to construct M trajectories of the solution (Y, Z) 
to the BSDE. 

In the following, the functions gi are chosen as step functions: 

T 

9i = l]t i - 1 ,t i ](t)/Vh, i = l,...,N, where h = — 

and the previous formulas are really simple (see (|10 p -(|li p and Proposition [S]) . Eventually, the 
main advantage of this method is that only one decomposition has to be computed per Picard 
iteration: the decomposition of F q . Therein lies the main difference between our approach and 
the approach based on regression technics developed by C. Bender and R. Denk in [BD07 . In 
their paper, for a given Picard iteration q and for each time U of the mesh grid, two projections 
have to be computed, one for Y t q and one for Z\. . The second difference comes from the way of 
computing Z q . In our method, once the decomposition of F q is computed, Z q is given explicitly 
as the Malliavin derivative of Y q . Let us also point out that our algorithm can handle fully path 
dependent terminal conditions. 

The rest of the paper is organized as follows. Section [2] contains the notations and the prelim- 
inary results, Section [3] describes precisely the algorithm, Section 0] is devoted to the study of the 
convergence of the algorithm and finally Section [5] contains some numerical experiments. Some 
technical proofs are post-done to the appendix. 



2 Preliminaries 

2.1 Definitions and Notations 

Given a probability space (f2, F, P) and an Revalued Brownian motion B, we consider 

• {{•^t)',t € [0,T]}, the filtration generated by the Brownian motion B and augmented 

• E t (X) denotes E(X\F t ) for any X in L 1 ^, F T , P). 

• L 2 (Ft) := L 2 (f2, Tti P) the space of all J-V-measurable random variables (r.v. in the follow- 
ing) X : n i — > R d satisfying ||X|| 2 = E(|X| 2 ) < oo. 

• S T (M. ) the space of all cadlag predictable processes cf> : fl x [0, T] i — > M d such that \\4>\\s2 = 
E ( su Pte[o.T] l^tl 2 ) < °°- 

• H T (R d ) the space of all predictable processes : Q x [0, T] i — > W 1 such that || 011^2 = 
E j^\<t>t\ 2 dt < co. 

• C^' 1 the set of continuously differentiable functions (j> : (t, x) £ [0, T] x R d with continuous 
and uniformly bounded derivatives w.r.t. t (resp. w.r.t. x) up to order k (resp. up to order 
I). The function (j> is also bounded. 

• ll^sp/lloo the norm of the derivatives of / : (t, x) £ [0, T] x M. w.r.t. all the space variables 
which sum equals j : Hd^/IlL : = E|fc|=j \\ d xl ' ' ' ^/IlL. whcrc 1^1 = h + ■ ■ ■ + k d . 

• the set of smooth functions / : ffi™ i — > R with partial derivatives of polynomial growth. 

• ||(-, -)!Il2 the norm on the space S T (M) x H T {R d ) defined by 

\\(Y,Z)\\l 2 =E[ sup (|F t | 2 )+ f T \Z t fdt). (3) 
\te[o,T] Jo I 
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We also recall some useful definitions related to Malliavin calculus. We use the notations of 
|Nua98j . 



• S denotes the class of random variables of the form F = f(W(hi),--- 1 W(h n )) 1 where 
v 



(h lr -- ,MeL 2 ([0,T];l") and W(h) = ^ h{t)dW t . 
• ll-Fjlr 2 denotes the following norm on S 



\F\\l 2 



:=E\F\ 2 + j2 E E (7 

«=l|a|=g \ J0 



1*1,— -tq) 



2 

dt\ ■ ■ ■ dt q 



where D a represents the multi-index Malliavin derivative operator. 
D r - 2 denotes the closure of S w.r.t. || • |L 2 and D 00 - 2 = n^ =1 W' 2 . 



2.2 Wiener Chaos Expansion 
2.2.1 Notations and useful results 

We refer to |Nua98j for more details on this section. Let us briefly recall the Wiener chaos expansion 
in the simple case of a real-valued Brownian motion. It is well known that every random variable 
F G L 2 (.Ft) as an expansion of the following form: 

F=E[F}+ f Ul ( Sl )dB Sl (4) 

u 2 (s 2 ,s 1 )dB Sl dB S2 + . . . + / •■•/ u n (s n ,...,si)dB sl ...dB Sn +... 

Jo Jo Jo 



where the functions (u n , n > 1) are deterministic functions. There is an ambiguity for the definition 
of these functions u„ . We adopt in this paper the following point of view: the function u n is defined 
on the simplex 

S n (T) := {( Sl , • • • , an) G [0, T] n : < Sl < . . . < s n < T} 
We define the iterated integral for a deterministic function / G L 2 (6>„(T)) as 

J n (f)~ f '/(sn,--- , Sl )dB Sl ---dB Sn . 

Jo Jo Jo 



10 Jo 

Due to the Ito isometry, ||J„(/)|| 2 = WfWh^^ and E[J n (f)J m (g)} = 5 nm < f,g >l=(s„(t))- 
Then, ||F|| 2 = E„>o H u «il^(5„(T))- 



Definition. Let F be a random variable in L 2 (.Ft) whose chaos expansion is given by ([!]). We 
introduce 

• P n {F) := J n (u n ) the Wiener chaos of order n of F. 

• C P (F) := Y^ n <pF'n(F) the chaos decomposition of F up to order p, i.e. 



C P (F) = E[F] + f Ul ( Sl )dB Sl + f [ 2 u 2 (s 2l s 1 )dB Sl dB s 
Jo Jo Jo 

r T r s P r s 2 
+ •••+/ / •••/ u p (s p ,...,s 1 )dB Sl ...dB Sp . (5) 
Jo "'O Jo 



We state two Lemmas useful for the sequel. 
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Lemma 1 (Nualart). F e B m ' 2 if and only if J2n>o( n + m- 1) X • • ■ x n x E[|P„(P)| 2 ] < oo. In 
this case, we have 

5^(n + m- 1) x ■ • • x n x E[|P„(P)| 2 ] < \\F\\l^ 2 . 

n>0 

From Lemma [TJ we deduce 
Lemma 2. Let F G D m > 2 . We have 

E[\F-C P (F)\ 2 } < 



(p + to) • • • (p + 1) 
Proo/. 

E[|P-C P (P)| 2 ] = £ E[P fc (P) 2 ] = ^ ( fc + m -l)..-fcx — — - xE[[P fc (P)p 

k>p+l k>p+l ^ ' 



□ 



The following Lemma gives some useful properties of the chaos decomposition. 
Lemma 3. Let F be a r.v. in L 2 (Pt) and H be in H|.(M). Then 
. Vp>l, E(|C P (P)| 2 ) <E(|P| 2 ) ; 

• c p(Jo H s ds ) = Jo C p (H s )ds. 

• For allt<r D t E r [C p (F)] = E r [C p -i(D t F)]. 

2.2.2 Wiener chaos expansion and Hermite polynomials 

Another approach to Wiener chaos expansion uses Hermite polynomials. This approach can be 
easily generalized when considering d-dimcnsional Brownian motions, this is then the one we 
consider in the following. We present it for d = 1. Let {gi}i>i be an orthonormal basis of 
L 2 (0,T). The Wiener chaos of order n, V n (F), is the L 2 -closure of the vector field spanned by 



\fn£.K ni ( J gj{s)dB s \ : \{nj)i>i\ ■= i 



i>l 

where K n is the Hermite polynomial of order n defined by the expansion: 

n>0 

with the convention K—\ = 0. With this normalization, we have K' n (x) = K n —i(x) for any integer 
n. It is well-known that (K n ) n >o is a sequence of orthogonal polynomials in L 2 (K, fx), where \i 
denotes the reduced centered Gaussian measure. Moreover, we have 

Kl(x)fi(dx) = — . 

Every square integrable random variable P, measurable with respect to Ft, admits the follow- 
ing orthogonal decomposition 



= (jf . 



(6) 
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where n = (rij)i>i is a sequence of positive integers and where \n\ stands for J^i>i n «- Taking into 
account the normalization of the Hermite polynomials we use, we get 

do=E[F], dl = n\ E F x JJ. >i A" rii N . 9l (s)dB^ , 

where n! = Y\i>i n i- Before describing the chaos decomposition formulas we use in the algorithm, 
we give a Lemma useful in the sequel. 

Lemma 4. Let g e L 2 (0, T) and let U t = J Q g 2 (s)ds. For n G N, let us define 
M? = U?' 2 K n (B(g) t /y/Ut) , B(g) t = J g{s)dB s . 

Then {M"}o<t<T is a martingale and 

dM? =g{t) M t n_1 dB t . 

2.3 Chaos decomposition formulas 

These formulas are based on the decomposition ©. To get tractable formulas, we consider a 
finite number of chaos and a finite number of functions (gx, ■ ■ ■ , giv)- The (gi)i<i<N functions 
are chosen such that we can quickly compute E(F\Ft) and D t E(F\FT) (as required in (0)). We 
develop in this Section the case d = 1, we refer to Section IB . 21 when d > 1. 



The first step consists in considering a finite number of chaos. In order to approximate the 
random variable F , we consider its projection C P (F) onto the first p chaos, namely 



c P (F) = *+E 1 < fe < p E |n|= ^ EL*- (j\(s)dB. 



(7) 



Of course, we still have an infinite number of terms in the previous sum and the second step 
consists in working with only the first N functions gi,. . . , g^ of an orthonormal basis of L 2 (0, T). 

Let us consider a regular mesh grid of N time steps T = {U = ijr,i = 0, • • • , TV} and the N 
step functions 

T 

9i = Ifa-uu] (t)/Vh, i = 1, ... ,N, where h := — . (8) 

We complete these N functions <?i,..., gn into an orthonormal basis of L 2 (0,T), ((?;);>i- For 
instance, one can consider the Haar basis on each interval (ti—i,U), i = 1,.. . , N. We implicitly 
assume that N > p. This leads to the following approximation 



c»{F) = do + E lnHk ^ ^ (/ T «(- . 



(9) 



where n = (ni, . . . , njv) and \n\ = m + . . . + nw- Due to the simplicity of the functions gj, 
i = 1, • • • , TV, we can compute explicitly 



/ gi{s)dB s = where d 

Jo 



B u - B ti 



Roughly speaking this means that Pk, the k th chaos, is generated by 

{K ni (Gi) . . . K nN (G N ) : m + . . . + n N = k} 
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Thus, the approximation we will use for the random variable F is 

C?(F) = d +E E d k R m (Gi) • • . Kn N (G N ) =do+J2J2<% II K ^ ( G *)> ( 10 ) 

fc=l|ri|=fc fc=l|ri|=fc l<i<N 

where the coefficients do and d£ are given by 

d = E[F], d% = n\E[FK ni (G 1 )...K nN (G N )\. (11) 

From (jllip . we deduce the expressions of E t (Cp F) and D t E t (C^(F)), useful for the approxi- 
mation of (Y,Z) by the chaos decomposition (see @). 

Proposition 5. Lei F be a real random variable in Ij 2 (J-t) and let r be an integer in {1, • ■ • , TV}. 

For all t r -i < t < t r , we have 

/ t — t r -\\ 2 / Bf — Bt r _ 1 



DA (<$-<*)) - *- V £ VJ 4 n (<r A',„( G ,, x (t^zl) Kn ,_, (*p*=t) , 

fc=l |n(r)|=fc V 7 V V 7 

n r >0 

where, if r < N and n = (ni, . . . , rijv), n(r) stands for (m, . . . , n r ). 

The proof of Proposition [5] is postponed to Section fB. II 
Remark 6. For t = t r , Proposition [5] leads to 

E tr (c»F) = d +j2 E ^LU^(^) 

fe=l | n (r)|=fc 
P 

D tr E ir (C*F) = /I" 1 / 2 ^ £ dl n j<r ^n, (Gi) x K nr -i (G r ) . 

fc=l | n (r)|=fc 

Tl r >0 

Let us end this subsection by some examples. 
Example 7 (Case p = 2). From ([TD |) -([TT |) . we have 

TV JV 3-1 AT 

Cf(F) - do + E^i(G,) + EE^^ 1 (G l )A' 1 (G J ) + E^ 2 (G,), 

3=1 j=l i=l 3=1 

where denotes the unit vector whose j component is one, and e%j = e% + ej. For j = 1, • ■ ■ , iV 
and z = 1 , • • • , j — 1 , it holds 

d? = E(FA' 1 (G J -)), = E(FA' 1 (G i )Jif 1 (G i )), d^ = 2E( J FA 2 (G J )). 
Remark [5] leads to 

7* r j—1 r 

E tr (CfF) = d + E rf i JA "i( G J ) +EE d 2 1 ^i(G,:)A- 1 (G,) +^d^A- 2 (G,), 

3=1 3=1 »=1 3=1 

D tr E tr (CfF) = /i- 1 / 2 K" +d 2e "A' 1 (G I .) +^d|^A'i(G i ) j . 
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3 Description of the algorithm 



The algorithm is based on three types of approximations : Picard's iterations, a Wiener chaos 
expansion up to a finite order and the truncation of a L 2 ([0, T}) basis in order to apply formulas 
of Proposition O We present the different steps of the approximation procedure in Section [3TTJ 
The practical implementation is presented in Section 13.21 



3.1 Approximation procedure 
3.1.1 Picard's iterations 

The first step consists in approximating (Y, Z) solution to (TTJ) by Picard's sequence 

(Y q , Z q ) q , built as follows : (Y° = 0, Z° = 0) and for all q > 1 



Y t q+1 = £ + j f (s, Y q , Zf) ds - £ Z q+1 • dB s , 0<t<T. 



(12) 



From ([13]), under the assumptions that £ G D 1 ' 2 and / G C°' 14 , we express (Y q+1 ,Z q+1 ) as a 
function of the processes (Y q , Z q ): 



7 +1 = Et U + £ f (*, Y q , Z q ) dsj , Z\ 



Y q+1 =l Z? +1 =D t Y q+1 , 

which can also be written 

Y t q+1 = Et U + jT f (s, Y q , Z q ) dsj J* f {a, Y q , Z q ) ds, Z q+1 = D t Y? +1 . (14) 

As recalled in the introduction, the computation of the conditional expectation is the corner- 
stone in the numerical resolution of BSDEs. Chaos decomposition formulas enable to circumvent 
this problem. 



3.1.2 Wiener Chaos Expansion 

Computing the chaos decomposition of the r.v. F = £ + f (s,Y q , Z q ) ds (appearing in (fT3")Q 
in order to compute Y t q+1 is not judicious. F depends on t, and then the computation of Y q+1 
on the grid T = {U = = 0, ■ • • , A^} would require N + 1 calls to the chaos decomposition 

function. To build a efficient algorithm, we need to call the chaos decomposition function as less 
as possible, since each call is computationally demanding and brings an approximation error due 
to the truncation and to the Monte-Carlo approximation (see next Sections). Then, we look for 
a r.v. F q independent of t such that Y t q+1 and Z q+1 can be expressed as functions of E t (F 9 ), 
D t E t (F q ) and of Y q and Z q . Equation (JT3J) gives a more tractable expression of Y q+1 . Let F q be 
defined by F q := £ + ^ f(s, Y q , Z q )ds. Then 

Y q+1 =E t (F q ) - f f( s ,Y q ,Z q )ds, Z q+1 = D t E t (F q ). (15) 
Jo 

The second type of approximation consists in computing the chaos decomposition of F q up to 
order p. Since F q does not depend on t, the chaos decomposition function C p is called only once 
per Picard's iteration. 

Let (Y q ' p , Z q ' p ) denote the approximation of (Y q , Z q ) built at step q using a chaos decompo- 
sition with order p: (Y *, Z ^) = (0, 0) and 

Y q+1 > p = E t [C p (F q 'P)} - f f (s, Y q 'P, ZfP) ds, Z q+1 > p = D t E t [C p {F™)\ , (16) 

Jo 

where F™ = £ + /„ / (s, Y™, Z™) ds. In the sequel, we also use the following equality 

Z q+1 ' p =E t [D t C p (F^)}. (17) 
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3.1.3 Truncation of the basis 

The third type of approximation comes from the truncation of the orthonormal L 2 ([0,T]) ba- 
sis used in the definition of C p (J7|)- Instead of considering a basis of L 2 ([0,T]), we only keep 
the first N functions (<?i, - • • , ffiv) defined by (J5J) to build the chaos decomposition function C p 
© • Proposition O gives us explicit formulas for E t (C™ F) and D t E t (C^F). Fr om p6p . we build 
((Y q ' p ' N , Z q ^ N ) q in the following way : ((Y°' P ' N , Z°^ N ) = (0,0) and 

yj«+i,p.# _ E t (Cp (F q ' p,N )) — [ f(s,Y q > p < N ,Z q ' p > N )ds, zf +x ' p ' N = D t (M t (Cp (F q,p ' N ))), (18) 

Jo 

where F q ' p < N := £ + J Q T f(s, Y q < p > N , Zf p > N )ds. 

Equation (fT5| is tractable as soon as we know closed formulas for the coefficients d k of the 
chaos decomposition oiE t (C* (F q ' p < N )) and D t (E t (Cp (F q ' p ' N ))) (see Proposition EJ). When it is 
not the case, we need to use a Monte-Carlo method to approximate these coefficients. The next 
Section is devoted to the practical implementation. In particular, we give the pseudo-code of the 
algorithm. 

3.2 Implementation 

In this Section, we first explain how to practically compute the chaos decomposition Cp (F) of a 
r.v. F . Then, we give the pseudo-code of the algorithm. 

3.2.1 Monte-Carlo simulations of the chaos decomposition 

Let F denote a r.v. of L 2 (JV). Practically, when we are not able to compute exactly do and/or 
the coefficients dZ of the chaos decomposition (|10p - (|ll|) of F, we use Monte-Carlo simulations to 
approximate them. Let (-F m )i< m <M be a M i.i.d. sample of F and (G™,- • ■ , G^)i<m<M be a 
M i.i.d. sample of (Gi,- • • ,Gn). We recall do and the coefficients (d%)i<k< p ,\ n \=k are given by 
d a = E[F] and d r k L = n\E[FK ni (Gi) . . . K nN (G N )] (see ((II])). Then, they are solutions of 

argmin E[\F - ip(c, G)| 2 ], (19) 

c = ( c >( c fc )l<fc<p ,|»| = fc) 

where tp : (c, G) i — > c + J2k=i J2\ n \=k c k Hi<i<n K n t (Gi). We propose two methods to approx- 
imate d := (d , (dfc)i<fc< P) | n |=jfe) 

• the first one consists in approximating the expectations of (1111) by empirical means (1m := 
( d "o,4i<fc< P ,|„| =fc ) where 

M . M 

* := M E pm > d l Yl S FmR nAG?) ■ --K nN (G%), 

m— 1 m— 1 

• the second one is based on a sample average approximation 

M 

d M := (do,d% 1<k<pH=k ) = argmin - £ |F m - ^(c, G™)| 2 

co,(c")i<fc<p,|„| = fc m=1 

Remark 8. In terms of computation time, the first method is much faster than the second one. 

• The first method requires 0(M x p) computations per coefficient. Since we are looking for 
0(N P ) coefficients, its computational cost is 0(M x p x N p ). 
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• The second method requires 0(M x p x N p ) computations to evaluate ^ Sm=i 1^™ — 
■0(c, G m )| 2 (in fact, it requires the same number of computations as the first method, since 
the function if; contains as much as additions as coefficients, and each addition contains as 
much as products as the associated coefficient). We still have to compute the argmin, which 
computational cost depends on the method we use. 

From a theoretical point of view, the second method gives better convergence results than 
the first one. For the first method, we only know that djyi converges to d a.s.. Concerning the 
second method, we know that dM converges to d a.s. and under regularity assumptions on ip, the 
uniform strong law of large numbers gives the a.s. convergence of -h Y^ m =i l-^™ ~~ ^(^Mj G m )\ 2 
to E[\F -ip{d,G)\ 2 }. 

In the following, Cp' M (F) denotes the approximation of the chaos decomposition of order p of 
F when using the first method to approximate the coefficients <ij!: 

C^ M (F) = do + J2 E ^ II K rn(Gi)- 

k=l \n\=k l<i<N 

E t (C^' M (F)) and D t {E t {Cp' M {F))) denote the conditional expectations obtained in Proposition 
El when (do,dfe)i<fc< P ,M=fc) are replaced by (do,d%)i<k< P ,\n\=k) '■ 

E, (C^F) :-i + ± « n, <r K„ (ft) * (tti) ^ A',. (^=f ) , 

D,MC; { F)):=H-V± Y. 5n i< ,A^(ft)x(^)^A V _ I (*=| if ), 

fc=l \n(r)\=k V 7 \ v r i / 

n r >0 

Remark 9. When M samples of C//' M (F) are needed, we can either use the same samples as the 

ones used to compute d and d£ : (C^(F)) m = d + Y7k=i J2\ n \=k d k Hi<i<N K m i G T)i or usc 
new ones. In the first case, we only require M samples of F and (Gi, • • • ,Gn)- The coefficients 
d£ and do are not independent of rii<i<7v ^"n, (G™)- The notation ~E t (Cp ,M (F)) introduced above 
cannot be linked to E (Cp' M F\Ft) ■ In the second case, the coefficients d% and do are independent 
of rii< l <Ar K nAG?) and we have E t (C£> M F) = E (C^ M F\F t ). This second approach requires 
2M samples of F and (Gi, • • • ,Gn) and its variance increases with TV. Practically, we use the 
first technique. 

3.2.2 Pseudo-code of the Algorithm 

In this Section, we describe in details the algorithm. We aim at computing M trajectories of an 
approximation of (Y, Z) on the grid T = {U = i = 0, • • • , N}. Starting from (Y°' P ' N , Z°'P' N ) = 
(0,0), (|18|) enables to get (Y q > p > N , Z q ' p ' N ) for each Picard's iteration q on T. However, if we 
only know the values of (Y q > p > N , Z q,p ' N ) on a grid and if we use a Monte Carlo procedure to 
compute the coefficients we are note able to compute (Y q+1 ' P ' N , Z q+1 ' P ' N ) on T exactly. Then, 
to take into account the different approximations presented before, we introduce F q ' p,N,M := 
^ + hE^o 1 m,Y q ^ M ,Z q f' N ' M ) and 

Y g+l,p,N,M = ^ ^N.M^NM^ - f (tj,Y^' P ' N ' M , Z q f> N ' M ) , 

3=0 

z q+i, P ,N,M = Du ( Eu ( C K<M(F q -P- N - M ))). (20) 



Here are the notations we use in the algorithm. 
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• d: dimension of the Brownian motion 

• q: index of Picard's iteration 

• Kn: number of Picard's iterations 

• M: number of Monte-Carlo samples 

• ./V: number of time steps used for the discretization of Y and Z 

• p: order of the chaos decomposition 

• Y q e Mn+i,mW represents M paths of Y q ' p ' N ' M computed on the grid T. 

• For all I G {1, • • ■ , d}, (Z c ')i e M n +im(^) represents M paths of (Z*<P' N ' M )i computed on 
the grid T . 

Since £ € L 2 (J-V), £ can be written as a measurable function of the Brownian path. Then, one 
gets one sample of £ from one sample of (G\, ■ • • , Gn) (where d represents B *' J^'^ 1 )■ 
For the sake of clearness, we detail the algorithm for d = 1. 



Algorithm 1 Iterative algorithm 

1: Pick at random N x M values of standard Gaussian r.v. stored in G. 

2: Using G, compute (£ m )i< m <M- 

3: Y° = 0, Z° = 0. 

4: for q = : K it - 1 do 
5: for m = 1 : M do 

6: Compute (F q ) m = f" + ^E,=i /(**, (Y q )i, m , (Z q ) t , m ) 

7: end for 

8: Compute the vector d = (do, (d% )i<k<p,\n\=k) °f the chaos decomposition of F q 
* d := i E"Li(F q ) m , d% = § YZli(F q ) m K n ,{G?) ■ ■ ■ K nN {G^) 

10: for j = : N - 1 do 

li: for m = 1 : M do 

12: Compute (E tj (C^ M F«))™, (D tj (E tj (C»> M F*))) m 

13: - (E tj (C^F*))™ + ELl /(**, (^'km, 

14: (Z^) J , m = (D ti (E t] (C»< M F q ))r 

15: end for 

16: end for 

17: end for 



Let us now deal with the complexity of the algorithm : 
For each q: 

• the computation of the vector F q (loop lineO requires 0(M x N) computations, 

• the computation of the vector d (line [8]) requires 0(M x p x (N x d) p ) computations, (in 
dimension d we have 0((N xd) p ) coefficients, and the computation of each coefficient requires 
0(M x p) computations (see Remark [5])), 

• for each N and M (lines fTDlITT] ! 

- the computation of (E t . (C^ M F q )) m and of (D l t , (E t . (Cg> M F q )))f^ d (line[I2J> requires 
0(d x p x (N x d) p ) computations 

— the computation of (Y 9+1 )j,m (loop linc ll3[) requires O(N) computations and the com- 
putation of ((Z q+1 ) l j m )i<i<d requires 0(d) computations. 

The complexity of the algorithm is then 0(K it x M x p x (N x d) p+1 ). 
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4 Convergence results 

We aim at bounding the error between (Y, Z) - the solution of ([T]) — and (Y q,p,N , Z q ' p ' N ) 
defined by ([2"TJ]) . Before stating the main result of the paper, we introduce some hypotheses. 

Hypothesis 10. Let F denote a r.v. in D'"' 2 such that for all integer r < m the function 

(s r ,-- - ►E(£>£>..., ar .F) 

is Holder of order a F , i.e. 3fcf s.t. \E(D^ t ... , Sr F) - E(dQ... u F)\ < fcf (|si — t\ | aF H + \s r - 

t r | QF ). In the following, denotes sup r<m k^ ■ 

The following Hypothesis is the generalization of Hypothesis [TU1 to the case m = oo. 

Hypothesis 11. Let F denote a r.v. in D 00,2 such that for all integer r the function 

is Holder of order a F , i.e. 3&f s.t. \E(D^ t ... . Sr F) - E(d{[\.. t F)\ < fef (|si - h \ a " + ■ • • + \s r - 

tr\ aF )- 

Theorem 12. Let k be an integer s.t. k < p. Assume that f G (j®> p > p an d £ satisfies Hypothesis 
\10\ when m = p. We have 

II Or - y**», z - z— < § + ^1 + A 2iq , P ) (1 

where Aq is given in Section \4-l\ A\ is given in Proposition \lb\ and Ai is given in Proposition 
If f G (^O' 00 ' 00 an( i ^ satisfies Hvvothesis \lll we get 

lim lim lim || (Y - Y q ' p ' N , Z - Z 9,P ' N ) ||? a = 0. 

Proof. We split the error in 3 terms : 

1. Picard's iterations : £ q = ||(Y - Y q , Z - Z q )\\^ 2 , where (Y q , Z q ) is defined by JT2J, 

2. the truncation of the chaos decomposition : £ q > p = \\(Y q - Y q > p ,Z q - Z q ' p )\\l 2 , where 
(Y9>p,Z«>p) is defined by (0, 

3. the truncation of the L 2 ([0,T]) basis : £ q ' p ' N = \\{Y q ' p - Y q ' p ' N , Z q > p - Z q ^ N )\\l 2 , where 
(Y q ' p ' N , Z q ' p ' N ) is defined by ([15]), 

We have 

\\(Y - Y q ' PlN , Z - Z q - p - N )\\l 2 < 3(£ q + £ q - p + £ q ' p ' N ). 
It remains to combine (|2~Tj) . Proposition [1d1 and Proposition |2"01 to get the first result. □ 

4.1 Picard's iterations 

The first type of error has already been studied in |PP92| and EPQ97| , we only recall the main 
result. 

Hypothesis 13. We assume 



the generator f : KxT — >R is Lipschitz continuous: there exists a constant Lj such that 
for all t G M+, y,z el andp,q G R d 

\f(t, y, z) - f(t,p, q)\ < L f (\y - P \ + \z- q\) , 
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• E[|£| 2 +/ T |/M,0)| 2 d s ]<^. 

From |EPQ97[ Corollary 2.1], we know that under Hypothesis [TBI the sequence (Y q ,Z q ) q 
defined by JTJJ) converges to (Y, Z) dP x dt a.s. and in S|(K) x Hf,(R d ). Moreover, we have 

E«;=\\(Y-Y«,Z-Z«)\\l 2 <^, (21) 
where A depends on T, ||£|| 2 and on ||/(-, 0, 0)|| 2 2 

[0,T] 

For the following, we also need 

Lemma 14. (\EPQ91\ Proof of Proposition 5.3]) Let to G W. Assume that f G c°' m ' m and 
£ G D m ' 2 . For all q £N, (Y q ,Z q ) belongs to L 2 ([0, T],B'™' 2 x (D m - 2 ) d ). 

In |EPQ97| , the proof of Lemma is done for m = 1, but it can be easily generalized for any 
integer m > 1. 

4.2 Error due to the truncation of the chaos decomposition 

We assume that the integrals are computed exactly, as well as expectations. The error is only due 
to the truncation of the chaos decomposition C p introduced in © . 

Lemma 15. Assume that f G C°' m ' m and £ G B m < 2 . For all q G N, (Y~«' p , Z q ' p ) belongs to 
L 2 ([0,T],B m < 2 x (B m - 2 ) d ). 

Proof of Lemma [73] Assume that (Y q ' p ,Z q ' p ) belongs to L 2 ([0, T];B m ' 2 x (B m < 2 ) d ). and let us 
show that (Y q+1 ' p , Z q+1 ' p ) belongs to L 2 ([0, T]; B m < 2 x (B m ' 2 ) d ). Since F q ' p G L 2 (.Ft), C p {F q - p ) G 
B m < 2 (in fact, we have C p {F q < p ) G B°°> 2 ). Then, E t [C p (F 9 ' p )] G B m < 2 . We deduce from Jig) that 
e L 2 ([0,T];B m - 2 ). It remains to prove that Z q+1 > p G L 2 ([0, T]; (B m - 2 ) d ). From (T7]) and 
the Clark-Ocone formula, we get f* Z q+1 ' p d B s = E t [C p (F q ' p )} - E[C p {F q - p )}. The r.h.s. belongs 
to B m - 2 . We get the result by using [PP921 Lemma 2.3], which proves that if an Ito integral is 
differcntiable in the Malliavin sense, its integrand is so. □ 

Proposition 16. Let to G N*. Assume that f G C°' m ' m and £ G B m < 2 . We recall E q < p = 
\\(Y q -Y q 'P,Z q - Z q -P)\\l 2 . We get 

gq +i,P < C T , T + i) L 2 £q , P + (22) 

where K x {m) depends on sup^^H d* p f H^), \\Z\\®m,i , T and on sup q£N \\{Y q , Z«)|| 2 2([0T] . Dm>2x(]D)m . 2)£i) 
and C\ is a scalar. 

Since £ 0,p = 0, we deduce from (|2"2"]1 that£ q ' p < ^^'^ where Ai(q,m) := K\ (to) 1 . 

Then, (Y p ' q , Z p ' q ) converges to (Y q ,Z q ) whenp tends to oo in ||(-, -)I|l 2 ('•see ([3J /or the Definition 
of the norm). 

Remark 17. We deduce from Proposition ll6l that for all T and L/, we have lim q ^ 00 lim p ^ f00 £ q ' p = 
0. When C\T(T + 1)L 2 < 1, i.e. for T small enough, we also get lim p ^ 00 lim q ^ 00 £ q ' p = 0. 

Proof of Proposition 1 1 61 For the sake of clearness, we assume d = 1 and to = 1. In the following, 
one notes AY t q ' p := Y q ' p - Y t q , AZf p := Z qp - Z q t and Af q t p ~ f{t,Y q *,Z q t p ) - f{t,Y q ,Z q ). 
Firstly, we deal with E[sup 0<t < T |Ay t 9+1 ' p | 2 ]. From CE5J) and (dHJ) we get 



AY q+hp =E t [C p (F q - p ) — F q ] — f Af q - p ds, 

Jo 



--E t [C p (0-t}+Et 



C p (£ f(s, Y q ' p , Z qp )ds^ - £ f(s, Y q , Z q )ds - J* Aff p ds 
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We introduce ±C P ( f(s, Y q , Zf)ds^ in the second conditional expectation. This leads to 



AY q+1 < p =E t [C p (Z)-$+E t 



Af q > p ds, 











+ E t 


[L 









C p (f(s,Y q ,Z q ))~f(s,Y q ,Z q )ds 



where we have used the second property of Lemma [3] to rewrite the third term. 

From the previous equation, we bound E[sup 0<t<T |AY t 9+1 ' p | 2 ] by using Doob's inequality and 
the Lipschitz property of / 



E[ sup \AY q+1 ' p \ 2 } < 16E[|C P (0 - CI 2 ] + 16E 

0<t<T 



C p / Affda 



16T / E 



\C p (f(s,Y s «,Z«))-f(s,Y q ,Z q )\< 



ds + 8TL 2 f 



|AF/' P | 2 + \AZ q s p \ 2 ]ds. 



To bound the second expectation of the previous inequality, we use the first property of Lemma [3] 
and the Lispchitz property of /. Then, wc bring together this term with the last one to get 



E[ sup |Ay t 9+1 ^| 2 ]<16E[|C p (0-C| 2 ] + 16T / E [|C p (/(s, F/, Z%)) — f(s, Yjf, Zf)\ 



0<t<T 



ds 



24TL 2 E[|AY S 9 < P | 2 + \AZf p \ 2 ]ds. 



(23) 



Let us now upper bound E[/ Q \AZj +1 *\ 2 ds\. To do so, wc use the Ito isometry E[ / Q \AZj +1 > p \ 2 ds\ 

E[( J Q T AZ$ +1 >PdB s ) 2 ]. Using the Definitions ([15 )1 - ([17 )1 of Z q+1 and Z q+1 ' p and the Clark-Oconc 
Theorem leads to 

f AZ q+1 < p dB s = F q - E(F q ) - (C p (F q ^) - E{C p (F q > p ))), 
Jo 

= Y q+1 + £ f(s, Y?, Z q )ds - Y q+1 - (Y q+1 - p + £ f(s, Y q 'P, Z™)ds - Y« +1 A 
Rearranging this summation makes appear AYp~ ' p and (AYq +1,p ). Young's inequality gives 

\AZ q+1 > p \ 2 ds <1e[su P \AY q+1 - p \ 2 ] + 32TL 2 f C ¥,[\AY q < p \ 2 + \AZ q < p \ 2 ]ds. (24) 

2 0<t<T Jo 



Since E[\AY^ P \ 2 + \AZf p \ 2 ]ds < (T+l)£^ p , by combining and (JUD we obtain 



-E q+1 * < 16E[|C P (£)-£| 2 ] + 16T J E[\C p (f(s,Y q ,Z«)) - f(s,Y q ,Z q )\ 2 \ ds + 56T(T + l)L 2 £ qp . 
Since £ G D 1 - 2 , / G C?' 1 ' 1 and (Y~« , Z c >) G L 2 ([0, T], D 1 ' 2 x D 1 ' 2 ) (see Lemma HU), Lemma gives 



2 

spj Moo 



+1 , p< 32 | ggrjjgg/ ii 

-^~PYII^IIl 2 ("x[o,t]) + 

112T(T+ 



|Ay 9 ||L2( Ox [ 0T ]) + || D t Z q \\1 2 ^ x [ ^dt 



Since (D t Y q , D t Z q ) converges to (D t Y,D t Z) in L 2 ([0, T]; Hf,(]R) x H|(K)) (see |EPQ97[ Proof of 
Proposition 5.3]), ^ follows. □ 
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4.3 Error due to the truncation of the basis 

We are now interested in bounding the error between (Y q ' p , Z q - p ) defined by Q16p and (Y q ' p ' , Z q ' p ' N ) 
denned by (|T8"j) . 

Before giving an upper bound for the error, we measure the error between C p and C p for a r.v. 
in D m ' 2 , when m> p. 

Remark 18. Let m £ N*, £ satisfies Hypothesis [TU1 and / G QO,m,m^ rp nell) f or a rj integers p and 
q, 7g iP := / Q T /(s, F s 9 ' p , Z qp )ds satisfies Hypothesis [TU1 where ai q p = \ A and ifm' p depends on 
#1, supfc^Jlld^/ljoo), T and on sup 9 ,< g || (Y q '> p , Z q ' ' p )\\ L 2 {[0iT] . nm .2 x[nm ,2 )d} . 
We refer to Section I A. 1 1 for the proof of Remark IT51 

Lemma 19. Let F denote a r.v. in L 2 (J-t) satisfying Hvvothesis 1 1 0\ for any integer m> p. We 
have 

E(\(C» - C P )(F)\ 2 ) < «) 2 (T^^j^i 2 ^ < {K F p f ^y F T(l + T)e T , 

where is defined in Hvvothesis ] 1(A 

We refer to Section [A. 21 for the proof of the Lemma. 
Proposition 20. Let m be an integer s.t. m > p. Assume that f G Q,™,,m an( ^ ^ satisfies 
HypothesisUE We recall S q ' p ' N := \\(Y q ' p - Y q < p > N \Z qp - Z q ' p ' N )\\l 2 . We get 

gq +i, P ,N < C2T{T + 1)L j £q ,p,N + K 2 ^ p ) 00 5 (25) 
where C 2 is a scalar and K2(q,p) depends on K^, sup fc<p (||9f : p /|| 2 <; ), T and on 

su P g '<<j \\(Y V ' P , Z q ' P )|lL2([o,T],Df' 2 x(DP 1 2 ) d )- 

Since £°> P > N = 0, we deduce from Q25J) that £ q > p ' N < A 2 (q,p) ) 2q « a \ where A 2 (q,p) := 
K 2 (q,p)T(T + ^ t [C c^t+1)l^-\ ■ Then > (Y p > q > N , Z p > q > N ) converges to (Y q > p ,Z q ' p ) when N 
tends to 00 in ||(-, -)||l 2 (see §3§ for the Definition of the norm). 



Proof of Proposition^^ For the sake of clearness, we assume d = 1. In the following, one notes 
AY?' P ' N := Y q ^ N - Y™, AZf p ' N := Z q ^ N - Zf p and Af q < p ' N := f(t,Y q ' p < N ,Z q ' p ' N ) - 
f(t,Y qp ,Z q ' p ). Firstly, we deal with E[sup < t < T \AY t q+1 ' p ' N \ 2 ]. From JTB]) and {THJ) we get 

AY t 9+1 ' p ' N = -E t [C£(F q ' p > N )-C p (F q > p )}+ f Af q ' p ' N ds. 

Jo 

Following the same steps as in the proof of Proposition [1(3 one gets 



E[ sup \AY q+1 ' p ' N \ 2 } <16E[|C^(0 -C p (0| 2 ] + 16E 



0<t<T 



K-c P ) ( j T f{. 



+ MTL) J E[\AY q ' p ' N \ 2 + \AZ q < p > N \ 2 ]ds. (26) 

Let us now upper bound E[J T \AZ q+1,p ' N \ 2 ds\. Following the same steps as in the proof of 
Proposition [TBI one gets 

/ \AZ q+r ' p ' N \ 2 ds <-E[ sup |Ay t ,+1 ' p,JV | 2 ] + 32TL? / E[\AY q < p > N \ 2 + \AZf p ' N \ 2 ]ds. 

Jo 2 o<t<T Jo 

(27) 
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Combining (f2"6"|) and ([2"T]) we obtain 



£ i+i,v,n <i 6E [|(C^ - C p )(0\ 2 } + 16E 



56T(T+ l)L 2 f £ q ' p ' N . 



(c 1 ; -Cp)( / q /(.s,rr,zr)ds 



Since £ and J 9jP := J f(s,Y q,p , Z q ' p )ds satisfy Hypothesis [TO] (see Remark IT8"| . Lemma [TO] gives 

/ rp\ 2aj Al 

f g+i, P ,JV < i 6 ( _ ) T(T+ l)e T ((/i«) 2 + (Kp«- P ) 2 ) + l\2T{T+l)L)£ qpN , 



and (S3 follows. 



□ 



5 Numerical Examples 



The computations have been done on a PC INTEL Core 2 Duo P9600 2.53 GHz with 4Gb of 
RAM. 



5.1 Non linear driver and path-dependent terminal condition 

We consider the case d=l, f(t 7 y, z) = cos(y) and £ = sup 0<t<T B t . 



• Convergence in p. Table Q] represents the evolution of Y r 



q,p,N,M 




and Z, 



q,p,N,M 




w.r.t q 



(Picard's iteration index), when p = 2 and p = 3. We fix M = 10 5 and N = 20. The seed 
of the generator is also fixed. 



iterations 


1 


2 


3 


4 


5 


6 


CPU time 


p = 2 


1.656357 


1.017117 


1.237135 


1.186691 


1.195462 


1.194256 


14.06 


p = 3 


1.656357 


1.012091 


1.234398 


1.183544 


1.192367 


1.191173 


174.09 


Table 1: Evolution of Y q ' p ' N ' M w.r.t. Picard's iterations, M = 10 5 , N = 20 and CPU time 


iterations 


1 


2 


3 


4 


5 


6 


CPU time 


p = 2 


0.969128 


0.249148 


0.525273 


0.459326 


0.470069 


0.469117 


14.06 


p = 3 


0.969128 


0.242977 


0.523846 


0.455827 


0.466903 


0.465939 


174.09 



Table 2: Evolution of Z, 



q,p,N,M 




.r.t. Picard's iterations, M = 10 5 , N = 20 and CPU time 



One notes that the difference between the values of Y r 



q,2,N,M 







and Y r 



q,3,N,M 







(resp. Zl 



q,2.N,M 



and Zq' ' ' ) doesn't exceed 0.2% (resp. 0.6%). This is due to the fast convergence of the 
algorithm in p. The CPU time is 12 times higher when p = 3 than when p = 2. Then, the 
use of order 3 in the chaos decomposition is not necessary. In the following, we take p = 2. 

Convergence in M. Figure [1] illustrates the evolution of Yq' p ' N ' M and Zq' p ' ' w.r.t. q 
whenp = 2 and N = 20 for different values of M. The seed of the generator is random. When 
M equals 10 4 and 10 5 the algorithm stabilizes after very few iterations. When M = 10 3 , 
there is no convergence. 

Convergence in N. Figure [5] illustrates the evolution of Yq' v ' ' m and Zq' p ' ' w.r.t. q 
when p = 2 and M = 10 5 for different values of N. The seed of the generator is random. 
The algorithm converges even when N = 5, but Yq' p ' 5 ' M is quite below Yq' p ' °' . 
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Evolution of Y"-q_0 Evolution of Z-^q O 




Figure 1: Evolution of Yq V ' ' and Z^' p ' ' w.r.t. q and M when TV = 20, p — 2 - £ = 
sup Q <t< T B t , f{t,y,z) = cos(y). 
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5.2 Linear Driver - Financial Benchmark 

We consider the case of pricing and hedging a Discrete Down and Out Barrier Call option, i.e. 
f{t, y, z) = —ry and £ := (St™ ^) + l\/ne[o,JV]s t >L: where S represents the Black-Scholes diffusion 

S t = S e (r ^ CT2)t+ffM \ Vt G [0,T]. 

The option parameters are r = 0.01, a = 0.2, T = 1, K = 0.9, X = 0.85, S = 1 and iV = 20 
(iV is also the number of time discretizations of the chaos decomposition). 



Evolution of Y_() Evolution of tlolta_() 




7 S 9 lO 1 1 12 13 14 "7 8 9 lO 11 12 13 14 



Figure 3: Evolution of r o 3 ' PlJVlJW and <5 := w.r.t. Zop(M) when iV = 20, p = 2, q = 5 

-Discrete Down and Out Barrier Call option 

We can get a benchmark for Yq and Zq by using a variance reduction Monte Carlo method. 
For this set of parameters, the reference values are Yq = 0.134267 with a confidence interval 

7.9468e-05 and S = = 0.8327. We compare them with Y Q q,p ' ' and \ So when N = 20, 
p = 2, q = 5 (we choose the first value of q from which the algorithm has converged) for different 
values of M. Figure [3] represents the evolution of Y" 5 ' p ' ' and <5q' p ' ' w.r.t. log(M). One 
notices that for M = 10 6 the computed values are very close to the reference ones. 

5.3 Non linear driver in dimension 5 - Financial Benchmark 

We consider the pricing and hedging of a Put Basket option in dimension 5, i.e. £ = (K — 

Vi = I,-- - ,5 S\ = S^'-^^' 713 *. 

fi l (resp. a 1 ) represents the trend (resp. the volatility) of the i th asset. B = (B 1 , • • • , B 5 ) is a 5- 
dimcnsional Brownian motion such that (_B l , B^) t = ptli^tj +tli—j. We suppose that p G (— -j, 1), 
which ensures that the matrix C = (plj^j + li=j)i<i,j<5 is positive definite. We also assume 
that the borrowing rate R is higher than the bond one r. In such a case, pricing and hedging 
the Put Basket option is equivalent to solving a BSDE with terminal condition £ and with driver 
/ defined by f(t,y,z) = -ry - 6 ■ z + (R - r)(y - ELi^^)*)" > whcrc 6 ■= - rl ) 

(1 is the vector whose every component is one) and S is the matrix defined by = cr l _Ly (L 
denote the lower triangular matrix involved in the Cholcsky decomposition C = LL*). We refer 
to |EPQ97| [Example 1.1] for more details. Figure H] represents the evolution of Y" 5 ' p ' ' , the 

approximated price at time 0, and the relative error on 5l := ^ — - - the quantity of 

asset 1 to possess at time — w.r.t. log(M). We compare our results with the ones obtained 
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log<IVl) 



log(M) 



Figure 4: Evolution of Yg' p ' N ' M and 8 (1) w.r.t. log(M) when N = 20, p = 2, g = 5, d = 5 - 
Basket Put option with different interest and borrowing rates 

using the Algorithm proposed in |GL10| (cited here as reference values) . The CPU time needed to 
compute price and delta when M = 50000 and N = 20 is 161s. One notices that the convergence 
is very fast and quite accurate for M = 50000. 

A Technical results of Section 14.3 
A.l Proof of Remark [181 

We prove the result for m = 1, i.e. we show that if £ satisfies Hypothesis [TU] for m = 1 and 
/ G C"' 1 ' 1 , then $:t — > E[D t J Q T f(s,Y™,Z™)ds] is Holder of order a := \ A a c with a 

constant K depending on K\, ||0 sp /||oo, T and on sup ? ,< 9 \\(Y q - p , Z q ■ p )\\L^([o,T];m^x(m.^)- Let 
us first prove the following Lemma. 

Lemma 21. Assume that £ satisfies Hvvothesis 1 1 0\ for m = 1 and f G C^' ' . For s,t such 
that < s < t and \t — s\ < 1, we have 

A??:=E[ sup (D t Y^' p - D s Y q,p ) 2 } + f E[{D t Z™ - £ s Z^) 2 ]dr < Ci(g,p)(t - s) 1A2a «. 

t<r<T Jt 

w/iere Ci(q,p) depends on K\, T, ||d sp /||oo and on sup g ,< g ||(T 9 '- p , Z 9 '' p )||l 2([0 T] D1 . 2x(D1 , 2)d) . 

Proof of Lemma Wf\ Let s, t, r be such that < s < t < r and \t — s\ < 1. Let us introduce some 
notations: A ts Y r <^ := D t Y™-D B Y™, A ts Z™ := D t Z™-D s Z™ and f(Qf) := f(r, Y™, Zf p ). 
From JTBJ and Lemmal we get A>7' p = E I .[C J) _i(Ai r9 ' p )] - D t f{B^ x )du. Then, 

By using Doob's inequality, Lemma[3]and the previous equation, we bound E[sup t<J , <T (At s Y' I ?' p ) 2 ]: 
E[ sup (A ts Y r ^) 2 } <12E[\D t F^ - D S F^\ 2 } 

t<r<T 

+ 3T fnDtf(6l l ) DJiOt^du + 3(t s) f E[\D s f{et l )\ 2 ]du. 

Jt Js 
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Using the definition of F** (see ([15])). we get D t F q ' p - D s F q ' p = D t £ - D s £ + ' D t f{O q ^ 1 ) 
D s f{9 q ^ 1 )du — J* D s f(9 q l ~ 1 )du. Plugging this result in the previous inequality yields 



E[ sup (A ts Y r ^) 2 } <36E[\D t t-DS 2 } 

t<r<T 



39T E[|A/(flT 1 )--D-/(flr 1 )| 2 ]du + 39(t-s) / E[\D s f {Ot^du 



Since = d y f {Q^OtY*- 1 * + d z f(6 g l ~ 1 )D t Z q - 1 'P and d y f and 0,/ arc bounded, we 

obtain 

E[ sup (A ts Y™) 2 } <36E[|Ae-^e| 2 ]+78T||a,/||L f E[| A^" 1 '^ 2 ]dw 

t<r<T Jt 

+ 78T||^/||L rEllA^'-^Hdu + TSCa^pJCt-*), (28) 



where C 2 (g,p) := K/IIL sup < s < T UA^'l^ + ||3 2 /||L sup < s < T ||£>.Z«* ||^ . Let us now 
upper bound E[\A ts Z^\ 2 ]dr. Using $T7J) and the Clark-Oconc formula gives J Q T Z™dB r = 
C p {Fi- l >P) - EiCpiFi- 1 *)). Hence, we have jf Z^dB r = Cpfi™- 1 *) - E t (C p (F9-bf)) = y™ + 
St f( e t l ) du - Y t" P - Thcn > since s<t<r,wc get 

^ A ts Z**dB r = A ts Y™ - A ts Y t q ' p + fiPtflfiT 1 ) - D s f(6l- l ))du. 

Young's inequality gives 

[ T E[\A ts Z^\ 2 ]dr < \e[ sup \A ts Y™\ 2 ] + 32T C E\\D t f{Q*- x ) - D s f{9t l )\] 2 du. 

Jt 1 t<r<T Jt 

As above, we develop the Malliavin derivatives of and we use that d y f and d z f are 

bounded. We obtain 

/ T E[|A ts Z^| 2 ]dr<iE[ sup | A ts F«f] + &4T\\dyf\\lo f E[| A ts Y«~ ^\ 2 }du (29) 

Jt 1 t<r<T Jt 

+ 64T||^/|| 2 ^ T E[|A ts zr 1 ^| 2 ]d M 
Combining (|28|) and (|29|) and using the Hypothesis [TO] satisfied by £ yields 

^A?* < (78C 2 (g,p) + 36(Kf) 2 )(t - s) 1A2a * + 284T||9 sp /|| 2 A?; 1 ^. 

Since A 4 'f = 0, we get Lemma I2"T1 by induction. □ 

We are now able to prove that $ is Holder. Let s, t be such that < s < t and \t — s\ < 1, we 
have 

$(t) - $(.s) = j\[D t f{6l) - D s f{6l)]dr - j' E{D s f(6«)]dr. 
As above, we develop the Malliavin derivatives of f{0 q ) 

$(t) - $(s) = j\[d y f{ei)A ts Y r ^]dr + j\[d z f(e q )A ts Z^}dr 
- f E[d y f{e*)D s Y r ™ + d z f(6?)D a Z™]dr. 
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Cauchy-Schwarz inequality gives E[d v f(0«)A ts Y r q > p } < ||9„/|| ooV /E[(A t8 r r ?,p ) 2 ]. Using the same 
argument to bound E[d z f(6*)At s Z?< p ] leads to 

Mi) - $( S )| 2 < 6T\\d sp f\\l,A™ + 6C 2 (q,p)(t - s), (30) 

where Af P (resp. C2(q,p)) has been introduced in Lemma |2"T1 (resp. in the proof of Lemma l2"Tj) . 
Combining ([50]) and Lemma |2"T1 ends the proof. 



A. 2 Proof of Lemma 1191 

We prove the result by induction. Lemma [TO1 is true for p = 0, since C^(F) = Cq(F). Assume 
that EdCC^ - C p ^)(F)n < {K^f (^f aF ^%. Since we have 

(C p N - C P )(F) = (C^ C P ^)(F) + (P p N - P P )(F), 
it remains to show that E(|(P p JV - P P )(F)\ 2 ) < (k£) 2 (f ) 2 " F p 2 |f . We recall 



P P (F) = [ [ ■ ■ ■ [ u p (s pi • • • , s 1 )dB Sl ■ ■ ■ dB Spl where u p : s p , ■ ■ ■ , s x i — ► E(L»^ ) .. D F), 
Jo Jo Jo Sp 

(31) 



P p( F )=J2 d % IT where d£ = n!E F ]J K tH (G t ) \ . (32) 

|n|=p l<i<AT y l<i<N J 

Let us rewrite P p (F) as a sum of stochastic integrals. Let r e N. Applying Lemma 0] to 
g : t i — >■ l ]ti _ liti] (t) yields M t r := h r / 2 A r ( ^"j!^" 1 ) is a martingale and M t r = M^~ 1 dB s . 
Then, M t r = J t *_ ■ • ■ / t ** M° ± dB Sl ■ ■ ■ dB Sr . For r = n, and t = < 4 , we get 



K ni (Gi 



1 



For \n\ := n\ 



K m {Gi) 



h 2 Jt t _ ± Jt 
n-N — P, we obtain 
1 



dB Sl ...dB a 



ti-i 



Ki<N 



d n = n' — 
d p n. hi 




dB 31 ---dB Sp , (33) 



n N integrals 



integrals m integrals 



N-l J tN-\ 



,(JV_l)| + 2 



*2 /"'infl 



ti A 



u p (l p , ■ ■ ■ ,h)dh ■ ■ - dip. (34) 



n N integrals n 2 integrals m integrals 

To compare P p (F) and P p N (F), we split the integrals in (|3"Tj) 



= E u 



s |n(N-l)|+2 /-*2 

ijv— i *i 



1(1)1+2 /Tl 

./0 



Up(s p , • ■ • , S!)dB 31 ■ ■■dB Sp . 



integrals 



n 2 integrals m integrals 



Combining 

l-T 



E 



and (123 yields E(|(P p JV - P P )(F)| 2 ) = 

S|n(JV-l)|+2 fta /•S|„(l)| +2 /•*! 



|n|= P r*^-J__^*"_i 



p Up\Sp, 

/l2 



Si) 



njy integrals 



integrals m integrals 



(35) 

i 

• • • dsp, 
(36) 
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Moreover, — § — u p (s p , ■ ■ • , si) 

h P Jt N -l Jt N -! Jt! Jt ± JO 



{u p (l p , ■ ■ ■ ,h) - Up(s p , • • • , si))dh ■ ■ - dl 



V 



n N integrals n 2 integrals m integrals 

Since u p satisfies Hypothesis [TU1 we get \u p (l p , ■ ■ ■ ,h) — u p (s p , ■ ■ ■ ,si)\ < k p {\l p — s p \ aF + • • • + 
|ii -si|" F ) < pkph a ". Then ^| - U p (s p , ■ ■ ■ , Si) <pk F h a ". Plugging this result in §jo$ ends 



the proof. 



B Wiener chaos expansion formulas 
B.l Proof of Proposition [5] 

Firstly, we compute E t (C^(F)) for t e]t r -i_,t r ]. From (JTHJ), we get 

p , 
E t (C»F) = d + Y,Y, d k Il <r K m(Gi) x Et (H> r K ^ G ^ 

k=l \n\=k 

Since Brownian increments are independent, we get Et r (Ili>r (Gi)) = ^n, (G r ) Yii> r ^[-^n; (^0] > 
which is null as soon as n r +i + • • • + tin > 0. Then, nested conditional expectations give 

p 

Et (CpF) = d a +J2 J2 d k H l<r K m ( G i) x E * ( R 'n r (Gr))- 

k=l \n(r)\=k 

By applying LemmaUwhen g : t \ — > l] tr _ 1>tr ](t), wegetE t (K nr (G r )) = f^-J K Ur f 
which yields the first result. Since K' n (x) = K n -±(x), the second result follows. 

B.2 Wiener chaos expansion formulas in M. d 

We want to approximate F € L 2 (J-t) using its chaos decomposition up to order p. We assume 
TV > dp. We consider the following truncated basis of L 2 ([0,T];]R d ) 

^U-UU]^) ■ , A7 • 1 , i , ^ 

= e, , z = 1, . . . , is, 7 = 1, .... a, where n= — 

where {t^ := i/i, i = 0, • • ■ , -/V} is a regular mesh grid and (&j)i<j<d represents the canonical basis 
of M. d . Pfe, the k th chaos, is generated by 



d N d N 



3 = 1 i=l o = l i=l 



For j = l,...,d, n J = (n\, . . . ,n J N ), one notes \n J \ = rv[ + . . . + n J N , n J \ = n\\ . . .n J N \ and for 

r < N, n j (r) = (n{, . . . ,n J r ). n = (n 1 , . . . , n d )* , \n\ = In 1 ! H h \n d \, n! = n 1 !...??.^ and 

n(r) = (n^r), . . . ,n d (r))*. Since the r.v. frii<j<d Ili<i<jv K rj. are orthogonal ones, the 

projection of F is given by 



k=l \n\=k 



i I ! 
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where the coefficients d£ are given by 

dfc = n\E IfTT TT A„ (g^ . 

L - l - L l<J<d-'--'-l<i<A r "i V v. 

Proposition 22. For t r _i < t < t r , we have 

fc=l |n(r)|=fc Krl<j<(! l<J<d V 7 \ V / 

and for I = 1, . . . , d, 

k=l \n{r)\=k 
n' T >0 

Remark 23. In particular, for t = t r and i = 1, . . . , d, 

fc=l | ra (r)|=fe 

A r (Et r (CpF)) = e e «*- i/a n t<r ru< d ^ ( g >< ^-i ( g ' ) k < H) ■ 

k=l |n(r)|=fe 
n l r >0 

Proof of Proposition We first compute E t (Cp F) for t £]t r _i,t r }. We have 

fc=l |n|=fc 

Since Brownian motions and their increments arc independents, wc get 

^ (iuru^ (<*)) = n.^Ai m n i>r n^ E 

which is null as soon as n\ +l + ■ ■ ■ + n l N + • • • + nf +1 + • • • + nfj > 0. Then, nested conditional 
expectations give 

k=l \n(r)\=k 

From Lemma HI for j = l,...,d M™ r := (t — tr-i)™*^ 2 K n i ^ ^ =i ^ is a martingale and 

dM"' = M"'" 1 l ]tr _ lttr] (t)dB{. Then, IIi<i<d (* - U-xf r ' 3 K n{ is also a martin- 

gale and the first result follows. Since K' , (x) = K n i _i(x), we get the second result. 

□ 

Conclusion. In this paper, we use Wiener chaos expansions together with the Picard pro- 
cedure to compute the solution to (H|). Once computed the chaos decomposition of F q , we get 
explicit formulas for both conditional expectations and the Malliavin derivative of conditional 
expectations. This enable to easily compute (Y q ,Z q ). Numerically, we obtain fast and accurate 
results, which encourage us to extend these results to other type of BSDEs, like 2-BSDEs. It is 
also possible to couple these Wiener chaos expansions together with the dynamic programming 
approach. This will be the subject of a forthcoming publication. 
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